--- title: "Large-Scale Drug Discovery: Enhanced QSAR Pipeline with 10,000+ Compounds" author: "Computational Drug Discovery - Production Pipeline" date: "`r Sys.Date()`" format: html: toc: true toc-depth: 3 code-fold: show code-tools: true theme: cosmo self-contained: true execute: warning: false message: false --- # 🧬 Production-Scale Drug Discovery Pipeline ## Project Overview \> 0.60 (vs. 0.21 in basic version)## Setup & Dependencies ```{r setup} #| label: setup #| message: false #Install Java 21 # In ~/.Rprofile or Rprofile.site Sys.setenv (JAVA_HOME = "C:/Program Files/Java/jdk-21" )Sys.setenv (PATH = paste (Sys.getenv ("PATH" ), "C:/Program Files/Java/jdk-21/bin/server" , sep = ";" ))# Core packages library (tidyverse)library (cli)library (glue)library (tictoc)library (progressr)library (future)library (furrr)# Cheminformatics library (rcdk)library (rJava)library (fingerprint)# Machine Learning library (ranger)library (xgboost)library (torch)library (luz)library (caret)# Visualization library (plotly)library (patchwork)library (ggridges)library (viridis)# Data retrieval library (httr)library (jsonlite)# Setup parallel processing plan (multisession, workers = parallel:: detectCores () - 1 )# Initialize Java for rcdk with maximum memory .jinit (parameters= "-Xmx8g" ) # Allocate 8GB RAM for Java :: cli_h1 ("🚀 Production-Scale Drug Discovery Pipeline Initialized" ):: cli_alert_info ("CPU Cores Available: {parallel::detectCores()}" ):: cli_alert_info ("GPU Available: {torch::cuda_is_available()}" )``` ## 1. Large-Scale Data Acquisition (10,000+ Compounds) ### 1.1 Enhanced Multi-Source Data Fetching ```{r data-acquisition-massive} #| label: fetch-massive-dataset :: cli_h2 ("Fetching Large-Scale Bioactivity Data from Multiple Sources" )# Enhanced function to fetch data with pagination and retry logic <- function (target_ids, max_compounds_per_target = 10000 ,activity_types = c ("IC50" , "Ki" , "Kd" , "EC50" )) {<- list ()for (target_id in target_ids) {:: cli_alert_info ("Processing target: {target_id}" )for (activity_type in activity_types) {:: cli_alert_info (" → Fetching {activity_type} data..." )<- 0 <- 1000 <- list ()repeat {# Rate limiting Sys.sleep (0.5 )<- glue ("https://www.ebi.ac.uk/chembl/api/data/activity.json?" ,"target_chembl_id={target_id}&" ,"standard_type={activity_type}&" ,"standard_relation==&" , # Exact measurements only "assay_type=B&" , # Binding assays "limit={limit}&" ,"offset={offset}" )<- tryCatch ({GET (url, timeout (60 ))error = function (e) {:: cli_alert_warning ("Request failed, retrying..." )Sys.sleep (5 )return (NULL )if (is.null (response) || status_code (response) != 200 ) {if (! is.null (response) && status_code (response) == 429 ) {:: cli_alert_warning ("Rate limit hit, waiting 60s..." )Sys.sleep (60 )next break <- tryCatch ({fromJSON (content (response, "text" , encoding = "UTF-8" ))error = function (e) {:: cli_alert_warning ("JSON parse error" )return (NULL )if (is.null (data) || length (data$ activities) == 0 ) break length (target_activities) + 1 ]] <- data$ activities %>% mutate (target_id = target_id, assay_type = activity_type):: cli_alert_info (" Retrieved {offset + nrow(data$activities)} compounds..." )<- offset + limitif (offset >= max_compounds_per_target || nrow (data$ activities) < limit) break if (length (target_activities) > 0 ) {<- bind_rows (target_activities)paste (target_id, activity_type, sep = "_" )]] <- combined:: cli_alert_success (" ✓ {nrow(combined)} {activity_type} records from {target_id}" )if (length (all_data) > 0 ) {<- bind_rows (all_data):: cli_alert_success ("Total raw records retrieved: {nrow(final_data)}" )return (final_data)else {return (NULL )# Comprehensive list of SARS-CoV-2 and related coronavirus targets <- c (# Primary SARS-CoV-2 targets "CHEMBL3927" , # Main protease (3CLpro) "CHEMBL5118" , # Mpro variant 1 "CHEMBL4523582" , # Mpro variant 2 "CHEMBL3927819" , # Related 3CLpro # Additional coronavirus proteases (for transfer learning) "CHEMBL2354684" , # SARS-CoV protease "CHEMBL4523580" , # MERS-CoV protease # Related viral proteases (structural similarity) "CHEMBL4523581" , # Additional coronavirus target "CHEMBL3927820" # Alternative 3CLpro isoform :: cli_alert_info ("Target compound goal: 10,000-20,000 unique molecules" ):: cli_alert_info ("Estimated fetch time: 30-60 minutes with rate limiting" )# Fetch comprehensive dataset <- fetch_chembl_comprehensive (target_ids = target_ids,max_compounds_per_target = 10000 ,activity_types = c ("IC50" , "Ki" ) # Focus on IC50 and Ki for consistency # Save raw data immediately dir.create ("data_massive" , showWarnings = FALSE )if (! is.null (bioactivity_raw_massive)) {<- bioactivity_raw_massive %>% select (where (~ is.atomic (.x) && ! is.matrix (.x)))write_csv (bioactivity_save, "data_massive/bioactivity_raw_massive.csv" )saveRDS (bioactivity_raw_massive, "data_massive/bioactivity_raw_massive.rds" ) # optional: keep original :: cli_alert_success ("Raw data saved: {nrow(bioactivity_save)} records" )``` ### 1.2 Advanced Data Preprocessing & Quality Control ```{r data-preprocessing-advanced} #| label: preprocess-massive-data :: cli_h2 ("Advanced Data Preprocessing & Quality Control" )<- function (bioactivity_raw) {:: cli_alert_info ("Starting preprocessing of {nrow(bioactivity_raw)} raw records..." )# Step 1: Flatten nested structures tryCatch ({<- bioactivity_raw %>% select (where (~ is.atomic (.x) && ! is.matrix (.x))):: cli_alert_success ("Step 1/6: Flattened nested columns ({ncol(bioactivity_flat)} columns)" )error = function (e) {:: cli_alert_danger ("Step 1 failed: {e$message}" )stop (e)# Step 2: Standardize units and calculate pActivity tryCatch ({<- bioactivity_flat %>% filter (! is.na (molecule_chembl_id),! is.na (canonical_smiles),! is.na (standard_value),! is.na (standard_units)) %>% mutate (standard_value_numeric = as.numeric (standard_value),# Convert all to nM standard_value_nm = case_when (== "nM" ~ standard_value_numeric,== "uM" ~ standard_value_numeric * 1000 ,== "pM" ~ standard_value_numeric / 1000 ,== "mM" ~ standard_value_numeric * 1e6 ,TRUE ~ NA_real_ pActivity = - log10 (standard_value_nm * 1e-9 ):: cli_alert_success ("Step 2/6: Standardized units to nM, calculated pActivity ({nrow(bioactivity_standardized)} records)" )error = function (e) {:: cli_alert_danger ("Step 2 failed: {e$message}" )print (head (bioactivity_flat))stop (e)# Step 3: Quality filters tryCatch ({<- bioactivity_standardized %>% filter (! is.na (pActivity),>= 4 ,<= 10 ,nchar (canonical_smiles) >= 5 ,nchar (canonical_smiles) <= 200 ,! grepl (" \\ ." , canonical_smiles) # Remove disconnected structures :: cli_alert_success ("Step 3/6: Applied quality filters (pActivity 4-10, valid SMILES) ({nrow(bioactivity_filtered)} records)" )error = function (e) {:: cli_alert_danger ("Step 3 failed: {e$message}" )stop (e)# Step 4: Handle duplicates intelligently tryCatch ({<- bioactivity_filtered %>% group_by (molecule_chembl_id, canonical_smiles) %>% summarise (pActivity = median (pActivity, na.rm = TRUE ),n_measurements = n (),targets = paste (unique (target_id), collapse = ";" ),assay_types = paste (unique (assay_type), collapse = ";" ),.groups = "drop" %>% filter (n_measurements >= 1 ):: cli_alert_success ("Step 4/6: Deduplicated (median pActivity per compound) ({nrow(bioactivity_dedup)} unique compounds)" )error = function (e) {:: cli_alert_danger ("Step 4 failed: {e$message}" ):: cli_alert_info ("Available columns: {paste(colnames(bioactivity_filtered), collapse = ', ')}" )stop (e)# Step 5: Activity classification tryCatch ({<- bioactivity_dedup %>% mutate (bioactivity_class = case_when (>= 7 ~ "Highly Active" ,>= 6 ~ "Active" ,TRUE ~ "Inactive" bioactivity_class = factor (bioactivity_class, levels = c ("Inactive" , "Active" , "Highly Active" )):: cli_alert_success ("Step 5/6: Classified activity levels ({nrow(bioactivity_classified)} records)" )# Debug: Print activity class counts <- table (bioactivity_classified$ bioactivity_class):: cli_alert_info ("Class counts: {paste(names(class_counts_debug), '=', class_counts_debug, collapse = ', ')}" )error = function (e) {:: cli_alert_danger ("Step 5 failed: {e$message}" )stop (e)# Step 6: Balance dataset (prevent class imbalance) tryCatch ({# Calculate max_per_class <- as.numeric (table (bioactivity_classified$ bioactivity_class))<- min (max (class_counts) * 1.5 , 5000 )<- as.integer (max_per_class):: cli_alert_info ("Max samples per class: {max_per_class}" )# Balance using group_map with .keep = TRUE to preserve grouping column <- bioactivity_classified %>% group_by (bioactivity_class) %>% group_map (~ {<- min (nrow (.x), max_per_class)%>% slice_sample (n = n_samples).keep = TRUE ) %>% # ADD THIS: .keep = TRUE preserves bioactivity_class bind_rows ():: cli_alert_success ("Step 6/6: Balanced classes (max {max_per_class} per class) ({nrow(bioactivity_balanced)} records)" )error = function (e) {:: cli_alert_danger ("Step 6 failed: {e$message}" ):: cli_alert_warning ("Skipping balancing step, returning unbalanced data" )<- bioactivity_classified# Final statistics tryCatch ({:: cli_h3 ("Final Dataset Statistics" ):: cli_alert_info ("Total unique compounds: {nrow(bioactivity_balanced)}" )# Verify bioactivity_class exists if (! "bioactivity_class" %in% colnames (bioactivity_balanced)) {:: cli_alert_danger ("bioactivity_class column missing!" ):: cli_alert_info ("Available columns: {paste(colnames(bioactivity_balanced), collapse = ', ')}" )return (bioactivity_balanced)<- bioactivity_balanced %>% group_by (bioactivity_class) %>% summarise (n = n (), .groups = "drop" ) %>% mutate (percentage = round (n / sum (n) * 100 , 1 )):: cli_alert_info ("Activity distribution:" )for (i in 1 : nrow (class_dist)) {:: cli_li ("{class_dist$bioactivity_class[i]}: {class_dist$n[i]} ({class_dist$percentage[i]}%)" )error = function (e) {:: cli_alert_warning ("Could not compute final statistics: {e$message}" )return (bioactivity_balanced)# Process the data with error handling :: cli_alert_info ("Starting preprocessing pipeline..." )<- tryCatch ({preprocess_massive_data (bioactivity_raw_massive)error = function (e) {:: cli_alert_danger ("Preprocessing failed completely: {e$message}" ):: cli_alert_info ("Traceback:" )print (traceback ())NULL # Save cleaned data if successful if (! is.null (bioactivity_clean_massive)) {tryCatch ({write_csv (bioactivity_clean_massive, "data_massive/bioactivity_clean_massive.csv" ):: cli_alert_success ("Cleaned data saved to data_massive/bioactivity_clean_massive.csv" ):: cli_alert_info ("Final dataset: {nrow(bioactivity_clean_massive)} compounds, {ncol(bioactivity_clean_massive)} columns" )error = function (e) {:: cli_alert_danger ("Failed to save CSV: {e$message}" ):: cli_alert_info ("Saving as RDS instead..." )saveRDS (bioactivity_clean_massive, "data_massive/bioactivity_clean_massive.rds" ):: cli_alert_success ("Saved as RDS file" )else {:: cli_alert_danger ("Preprocessing failed - no data to save" )``` ### 1.3 Data Visualization Dashboard ```{r data-viz-dashboard} #| label: data-visualization #| fig-width: 16 #| fig-height: 12 :: cli_h3 ("Generating Data Visualization Dashboard" )# 1. pActivity Distribution <- plot_ly (bioactivity_clean_massive) %>% add_histogram (x = ~ pActivity,marker = list (color = "#3498db" ,line = list (color = "white" , width = 1 )),nbinsx = 50 ) %>% layout (title = paste ("pActivity Distribution (n =" , nrow (bioactivity_clean_massive), "compounds)" ),xaxis = list (title = "pActivity" ),yaxis = list (title = "Count" ))# 2. Activity Class Distribution <- bioactivity_clean_massive %>% group_by (bioactivity_class) %>% summarise (n = n (), .groups = "drop" )<- plot_ly (class_counts) %>% add_pie (labels = ~ bioactivity_class,values = ~ n,marker = list (colors = c ("#e74c3c" , "#f39c12" , "#27ae60" )),textinfo = "label+value+percent" ,hole = 0.4 ) %>% layout (title = "Activity Class Distribution" )# 3. Density Ridge Plot <- ggplot (bioactivity_clean_massive, aes (x = pActivity, y = bioactivity_class, fill = bioactivity_class)) + geom_density_ridges (alpha = 0.7 , scale = 2 ) + scale_fill_manual (values = c ("#e74c3c" , "#f39c12" , "#27ae60" )) + theme_minimal () + labs (title = "pActivity Distribution by Class (Ridge Plot)" ,x = "pActivity" , y = "" ) + theme (legend.position = "none" ,text = element_text (size = 12 ))ggplotly (fig3)``` ## 2. Parallel Molecular Descriptor Computation ### 2.1 Parallel ECFP4 Fingerprints ```{r ecfp-parallel} #| label: ecfp-parallel-computation :: cli_h2 ("Computing ECFP4 Fingerprints (Parallel Processing)" )<- function (smiles_vector, radius = 2 , nbits = 1024 ) {:: cli_alert_info ("Computing ECFP{radius*2} for {length(smiles_vector)} molecules (parallel)..." )# Split into chunks for parallel processing <- parallel:: detectCores () - 1 <- ceiling (length (smiles_vector) / n_chunks)<- split (smiles_vector, ceiling (seq_along (smiles_vector) / chunk_size))# Parallel computation with progress handlers ("cli" )with_progress ({<- progressor (steps = length (chunks))<- future_map (chunks, function (chunk) {p ()map_dfr (chunk, function (smi) {tryCatch ({<- parse.smiles (smi)[[1 ]]<- get.fingerprint (mol, type = "circular" , fp.mode = "bit" ,circular.type = "ECFP4" , size = nbits)<- as.vector (fp@ bits)<- integer (nbits)<- 1 as.data.frame (t (fp_binary))error = function (e) {as.data.frame (t (integer (nbits))).options = furrr_options (seed = TRUE ))<- bind_rows (ecfp_chunks)colnames (ecfp_matrix) <- paste0 ("ECFP" , radius* 2 , "_" , 1 : nbits):: cli_alert_success ("ECFP{radius*2} computation complete!" )return (ecfp_matrix)tic ("ECFP4 Parallel Computation" )<- compute_ecfp_parallel (bioactivity_clean_massive$ canonical_smiles,radius = 2 , nbits = 1024 )toc ()``` ### 2.2 Parallel MACCS Keys ```{r maccs-parallel} #| label: maccs-parallel-computation :: cli_h2 ("Computing MACCS Keys (Parallel Processing)" )<- function (smiles_vector) {:: cli_alert_info ("Computing MACCS keys for {length(smiles_vector)} molecules (parallel)..." )<- parallel:: detectCores () - 1 <- ceiling (length (smiles_vector) / n_chunks)<- split (smiles_vector, ceiling (seq_along (smiles_vector) / chunk_size))handlers ("cli" )with_progress ({<- progressor (steps = length (chunks))<- future_map (chunks, function (chunk) {p ()map_dfr (chunk, function (smi) {tryCatch ({<- parse.smiles (smi)[[1 ]]<- get.fingerprint (mol, type = "maccs" )<- as.vector (fp@ bits)<- integer (166 )<- 1 as.data.frame (t (maccs_binary))error = function (e) {as.data.frame (t (integer (166 ))).options = furrr_options (seed = TRUE ))<- bind_rows (maccs_chunks)colnames (maccs_matrix) <- paste0 ("MACCS_" , 1 : 166 ):: cli_alert_success ("MACCS computation complete!" )return (maccs_matrix)tic ("MACCS Parallel Computation" )<- compute_maccs_parallel (bioactivity_clean_massive$ canonical_smiles)toc ()``` ### 2.3 Enhanced Lipinski Descriptors (Parallel) ```{r lipinski-parallel} #| label: lipinski-parallel-computation # ============================================================================ # 2.3 Enhanced Lipinski Descriptors (PROVEN SERIAL VERSION - WORKS!) # ============================================================================ :: cli_h2 ("Computing Enhanced Lipinski Descriptors (Serial - Reliable)" )<- function (smiles_vector) {:: cli_alert_info ("Computing Lipinski descriptors for {length(smiles_vector)} molecules..." )# Parse SMILES first :: cli_alert_info ("Step 1: Parsing SMILES structures..." )<- purrr:: map (smiles_vector, function (smi) {tryCatch ({parse.smiles (smi)[[1 ]]error = function (e) NULL )<- sum (! sapply (molecules, is.null)):: cli_alert_info ("Successfully parsed {successful_parses}/{length(molecules)} SMILES" )# Progress bar <- length (molecules)<- txtProgressBar (min = 0 , max = total, style = 3 , title = "Lipinski Progress" )# Calculate descriptors serially with better error handling <- list ()for (i in seq_along (molecules)) {setTxtProgressBar (pb, i)<- molecules[[i]]if (is.null (mol)) {<- data.frame (MW = NA_real_ , LogP = NA_real_ , HBD = NA_integer_ , HBA = NA_integer_ ,TPSA = NA_real_ , nRotB = NA_integer_ , nAtoms = NA_integer_ ,Aromatic_Bonds = NA_integer_ , Ring_Count = NA_integer_ next # Try each descriptor individually with error catching tryCatch ({<- tryCatch (get.mol2formula (mol)@ mass, error = function (e) NA_real_ )<- tryCatch (eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.XLogPDescriptor" )[[1 ]], error = function (e) NA_real_ )<- tryCatch (eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.HBondDonorCountDescriptor" )[[1 ]], error = function (e) NA_integer_ )<- tryCatch (eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.HBondAcceptorCountDescriptor" )[[1 ]], error = function (e) NA_integer_ )<- tryCatch (eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.TPSADescriptor" )[[1 ]], error = function (e) NA_real_ )<- tryCatch (eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.RotatableBondsCountDescriptor" )[[1 ]], error = function (e) NA_integer_ )<- tryCatch (get.atom.count (mol), error = function (e) NA_integer_ )<- tryCatch (eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.AromaticBondsCountDescriptor" )[[1 ]], error = function (e) NA_integer_ )<- tryCatch (eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.RingCountDescriptor" )[[1 ]], error = function (e) NA_integer_ )<- data.frame (MW = mw, LogP = logp, HBD = hbd, HBA = hba,TPSA = tpsa, nRotB = nrotb, nAtoms = natoms,Aromatic_Bonds = aromatic_bonds, Ring_Count = ring_counterror = function (e) {<<- data.frame (MW = NA_real_ , LogP = NA_real_ , HBD = NA_integer_ , HBA = NA_integer_ ,TPSA = NA_real_ , nRotB = NA_integer_ , nAtoms = NA_integer_ ,Aromatic_Bonds = NA_integer_ , Ring_Count = NA_integer_ close (pb)# Combine all results <- bind_rows (lipinski_list)# Report statistics <- sapply (lipinski_data, function (x) sum (is.na (x))):: cli_alert_success ("Lipinski descriptor calculation complete!" ):: cli_h3 ("NA Summary per Descriptor:" )for (col in names (na_summary)) {:: cli_li ("{col}: {na_summary[col]} NAs ({round(100*na_summary[col]/nrow(lipinski_data), 1)}%)" )return (lipinski_data)# Execute with timing tic ("Lipinski Serial Computation" )<- compute_lipinski_serial_reliable (bioactivity_clean_massive$ canonical_smiles)toc ()# Save immediately write_csv (lipinski_massive_new, "data_massive/lipinski_descriptors_massive_fixed.csv" ):: cli_alert_success ("Fixed Lipinski descriptors saved!" )``` ### 2.4 Combine All Descriptors ```{r} # Check NA counts per descriptor type #cli::cli_alert_info("NA counts:") #cli::cli_alert_info("Lipinski NAs: {sum(is.na(lipinski_massive_new))}") #cli::cli_alert_info("MACCS NAs: {sum(is.na(maccs_fps_massive))}") #cli::cli_alert_info("ECFP4 NAs: {sum(is.na(ecfp4_fps_massive))}") # Check which molecules have all NA descriptors #lipinski_na_rows <- rowSums(is.na(lipinski_massive_new)) == ncol(lipinski_massive_new) #cli::cli_alert_info("Molecules with all NA Lipinski: {sum(lipinski_na_rows)}") ``` ```{r combine-descriptors-massive} #| label: combine-massive-features #cli::cli_h2("Combining All Molecular Descriptors") #full_dataset_massive <- bioactivity_clean_massive %>% # bind_cols(lipinski_massive_new) %>% # bind_cols(maccs_fps_massive) %>% # bind_cols(ecfp4_fps_massive) %>% # select(-targets, -assay_types) %>% # drop_na() #cli::cli_alert_success("Combined dataset: {nrow(full_dataset_massive)} compounds") #cli::cli_alert_success("Total features: {ncol(full_dataset_massive) - 4}") #write_csv(full_dataset_massive, "data_massive/full_dataset_massive.csv") ``` ```{r} :: cli_h2 ("Fixing Ring Count Descriptor (Post-Processing)" )# Load the Lipinski data you just computed <- read_csv ("data_massive/lipinski_descriptors_massive_fixed.csv" )# Fix Ring_Count by computing it a different way (simpler approach) <- function (smiles_vector) {:: cli_alert_info ("Computing Ring_Count safely for {length(smiles_vector)} molecules..." )<- txtProgressBar (min = 0 , max = length (smiles_vector), style = 3 )<- purrr:: map_dbl (seq_along (smiles_vector), function (i) {setTxtProgressBar (pb, i)<- smiles_vector[i]tryCatch ({<- parse.smiles (smi)[[1 ]]if (is.null (mol)) return (NA_real_ )# Try RDKit Ring descriptor eval.desc (mol, "org.openscience.cdk.qsar.descriptors.molecular.RingCountDescriptor" )[[1 ]]error = function (e) {NA_real_ # Return NA if fails close (pb)return (ring_counts)# Compute Ring_Count <- compute_ring_count_safe (bioactivity_clean_massive$ canonical_smiles)# Count how many we recovered <- sum (is.na (ring_counts)):: cli_alert_info ("Ring_Count computation: {length(ring_counts) - na_count}/{length(ring_counts)} successful ({round(100*(1-na_count/length(ring_counts)), 1)}%)" )# If Ring_Count is still mostly NA, just remove it or fill with zeros if (na_count > length (ring_counts) * 0.5 ) {:: cli_alert_warning ("Ring_Count has >50% NAs, removing this descriptor" )$ Ring_Count <- NULL else {$ Ring_Count <- ring_counts:: cli_alert_success ("Lipinski descriptors fixed!" )# Now combine EVERYTHING with the original names :: cli_h2 ("Combining All Molecular Descriptors" )<- bioactivity_clean_massive %>% bind_cols (lipinski_massive) %>% bind_cols (maccs_fps_massive) %>% bind_cols (ecfp4_fps_massive) %>% select (- targets, - assay_types) %>% drop_na ():: cli_alert_success ("✅ Combined dataset: {nrow(full_dataset_massive)} compounds" ):: cli_alert_success ("✅ Total features: {ncol(full_dataset_massive) - 4}" )# Feature breakdown <- ncol (lipinski_massive)<- 166 <- 1024 :: cli_alert_info ("Feature breakdown:" ):: cli_li ("Lipinski: {n_lipinski} descriptors" ):: cli_li ("MACCS: {n_maccs} keys" ):: cli_li ("ECFP4: {n_ecfp4} fingerprints" ):: cli_li ("TOTAL: {n_lipinski + n_maccs + n_ecfp4} features" )# Save with original name write_csv (full_dataset_massive, "data_massive/full_dataset_massive.csv" ):: cli_alert_success ("🎉 Full dataset saved!" )``` ## 3. Optimized Train/Test Split ```{r train-test-split-massive} #| label: data-split-massive :: cli_h2 ("Creating Stratified Train/Test Split (80/20)" )set.seed (42 )<- createDataPartition (full_dataset_massive$ pActivity, p = 0.8 , list = FALSE )<- full_dataset_massive[train_indices, ]<- full_dataset_massive[- train_indices, ]:: cli_alert_success ("Training set: {nrow(train_data_massive)} compounds" ):: cli_alert_success ("Test set: {nrow(test_data_massive)} compounds" )write_csv (train_data_massive, "data_massive/train_data_massive.csv" )write_csv (test_data_massive, "data_massive/test_data_massive.csv" )``` ## 4. Optimized Machine Learning Models ### 4.1 Random Forest (Production-Optimized) ```{r random-forest-production} #| label: rf-production :: cli_h2 ("Training Production Random Forest Model" )<- setdiff (names (train_data_massive),c ("molecule_chembl_id" , "canonical_smiles" ,"pActivity" , "bioactivity_class" , "n_measurements" )):: cli_alert_info ("Training on {length(feature_cols_massive)} features..." )tic ("Random Forest Production Training" )<- ranger (~ .,data = train_data_massive %>% select (pActivity, all_of (feature_cols_massive)),num.trees = 1000 , # Increased trees for large dataset mtry = floor (sqrt (length (feature_cols_massive))),min.node.size = 5 ,importance = "permutation" ,num.threads = parallel:: detectCores () - 1 ,verbose = TRUE ,seed = 42 toc ()# Predictions <- predict (rf_model_massive, test_data_massive)$ predictions<- data.frame (Actual = test_data_massive$ pActivity,Predicted = rf_predictions_massive,Residuals = test_data_massive$ pActivity - rf_predictions_massive<- sqrt (mean (rf_metrics_massive$ Residuals^ 2 ))<- cor (rf_metrics_massive$ Actual, rf_metrics_massive$ Predicted)^ 2 <- mean (abs (rf_metrics_massive$ Residuals)):: cli_alert_success ("Random Forest - RMSE: {round(rf_rmse, 4)}, R²: {round(rf_r2, 4)}, MAE: {round(rf_mae, 4)}" )dir.create ("results_massive" , showWarnings = FALSE )write_csv (rf_metrics_massive, "results_massive/rf_predictions.csv" )saveRDS (rf_model_massive, "results_massive/rf_model_massive.rds" )``` ### 4.2 XGBoost (GPU-Accelerated) ```{r xgboost-production} #| label: xgb-production # ============================================================================ # STEP 1: PREPARE DATA FOR XGBOOST # ============================================================================ :: cli_h2 ("Preparing XGBoost Data Matrices" )# Define feature columns (exclude metadata and target) <- setdiff (names (train_data_massive), c ("molecule_chembl_id" , "canonical_smiles" , "pActivity" , "n_measurements" )):: cli_alert_info ("Number of features: {length(feature_cols_massive)}" )# Create training data and CONVERT TO NUMERIC <- as.matrix (train_data_massive %>% select (all_of (feature_cols_massive)) %>% mutate (across (everything (), as.numeric)))<- as.numeric (train_data_massive$ pActivity)# Create test data and CONVERT TO NUMERIC <- as.matrix (test_data_massive %>% select (all_of (feature_cols_massive)) %>% mutate (across (everything (), as.numeric)))<- as.numeric (test_data_massive$ pActivity):: cli_alert_info ("Data types check:" ):: cli_li ("train_x_massive class: {class(train_x_massive)}" ):: cli_li ("train_x_massive mode: {mode(train_x_massive)}" ):: cli_li ("train_y_massive class: {class(train_y_massive)}" )# Convert to xgboost DMatrix format <- xgb.DMatrix (data = train_x_massive, label = train_y_massive)<- xgb.DMatrix (data = test_x_massive, label = test_y_massive):: cli_alert_success ("✅ XGBoost matrices created successfully!" ):: cli_li ("Training matrix: {nrow(train_x_massive)} rows × {ncol(train_x_massive)} features" ):: cli_li ("Test matrix: {nrow(test_x_massive)} rows × {ncol(test_x_massive)} features" )# ============================================================================ # STEP 2: TRAIN XGBOOST MODEL # ============================================================================ :: cli_h2 ("Training XGBoost Model (CPU)" )<- list (objective = "reg:squarederror" ,eval_metric = "rmse" ,max_depth = 8 ,eta = 0.05 ,subsample = 0.8 ,colsample_bytree = 0.8 ,min_child_weight = 5 ,gamma = 0.1 ,tree_method = "hist" ,nthread = parallel:: detectCores () - 1 tic ("XGBoost Training (CPU)" )<- xgb.train (params = params_massive,data = dtrain_massive,nrounds = 2000 ,watchlist = list (train = dtrain_massive, test = dtest_massive),early_stopping_rounds = 100 ,verbose = 1 ,print_every_n = 100 toc ():: cli_alert_success ("✅ XGBoost training complete!" )# ============================================================================ # STEP 3: GET PREDICTIONS & METRICS # ============================================================================ <- predict (xgb_model_massive, dtest_massive)<- data.frame (Actual = test_y_massive,Predicted = xgb_predictions_massive,Residuals = test_y_massive - xgb_predictions_massive<- sqrt (mean (xgb_metrics_massive$ Residuals^ 2 ))<- cor (xgb_metrics_massive$ Actual, xgb_metrics_massive$ Predicted)^ 2 <- mean (abs (xgb_metrics_massive$ Residuals)):: cli_alert_success ("XGBoost Results:" ):: cli_li ("RMSE: {round(xgb_rmse, 4)}" ):: cli_li ("R²: {round(xgb_r2, 4)}" ):: cli_li ("MAE: {round(xgb_mae, 4)}" )# ============================================================================ # STEP 4: SAVE RESULTS # ============================================================================ write_csv (xgb_metrics_massive, "results_massive/xgb_predictions.csv" )xgb.save (xgb_model_massive, "results_massive/xgb_model_massive.json" ):: cli_alert_success ("✅ XGBoost model and predictions saved!" )``` ### 4.3 Deep Neural Network (GPU-Optimized, Batch Training) ```{r dnn-production} #| label: dnn-production :: cli_h2 ("Training Production Deep Neural Network (GPU-Optimized)" )# Prepare data # Select feature columns <- train_data_massive %>% select (all_of (feature_cols_massive))# Convert all columns to numeric (handles factors, characters) <- train_features_df %>% mutate (across (everything (), as.numeric))# Convert to numeric matrix for model input <- as.matrix (train_features_num)# Target variable <- train_data_massive$ pActivity# For test set similarly <- test_data_massive %>% select (all_of (feature_cols_massive))<- test_features_df %>% mutate (across (everything (), as.numeric))<- as.matrix (test_features_num)<- test_data_massive$ pActivity# Normalize <- colMeans (train_x_massive)<- apply (train_x_massive, 2 , sd)== 0 ] <- 1 # Prevent division by zero <- scale (train_x_massive, center = mean_vals_massive, scale = sd_vals_massive)<- scale (test_x_massive, center = mean_vals_massive, scale = sd_vals_massive)# Convert to tensors <- if (torch:: cuda_is_available ()) "cuda" else "cpu" :: cli_alert_info ("Using device: {device}" )<- torch_tensor (train_x_norm_massive, dtype = torch_float32 ())$ to (device = device)<- torch_tensor (train_y_massive, dtype = torch_float32 ())$ unsqueeze (2 )$ to (device = device)<- torch_tensor (test_x_norm_massive, dtype = torch_float32 ())$ to (device = device)<- torch_tensor (test_y_massive, dtype = torch_float32 ())$ unsqueeze (2 )$ to (device = device)# Production DNN architecture <- nn_module ("TabularResNet" ,initialize = function (input_dim) {$ fc1 <- nn_linear (input_dim, 1024 )$ bn1 <- nn_batch_norm1d (1024 )$ drop1 <- nn_dropout (0.4 )$ fc2 <- nn_linear (1024 , 1024 )$ bn2 <- nn_batch_norm1d (1024 )$ drop2 <- nn_dropout (0.3 )$ fc3 <- nn_linear (1024 , 512 )$ bn3 <- nn_batch_norm1d (512 )$ drop3 <- nn_dropout (0.3 )$ fc4 <- nn_linear (512 , 256 )$ bn4 <- nn_batch_norm1d (256 )$ drop4 <- nn_dropout (0.2 )$ fc5 <- nn_linear (256 , 128 )$ bn5 <- nn_batch_norm1d (128 )$ drop5 <- nn_dropout (0.2 )$ fc6 <- nn_linear (128 , 64 )$ bn6 <- nn_batch_norm1d (64 )$ drop6 <- nn_dropout (0.1 )$ output <- nn_linear (64 , 1 )$ relu <- nn_relu ()forward = function (x) {<- x<- self$ fc1 (x0)<- self$ bn1 (x)<- self$ relu (x)<- self$ drop1 (x)# Residual block 1 <- x<- self$ fc2 (x)<- self$ bn2 (x)<- self$ relu (x)<- self$ drop2 (x)<- x + x_skip # Residual # Feedforward blocks <- self$ fc3 (x)<- self$ bn3 (x)<- self$ relu (x)<- self$ drop3 (x)<- self$ fc4 (x)<- self$ bn4 (x)<- self$ relu (x)<- self$ drop4 (x)<- self$ fc5 (x)<- self$ bn5 (x)<- self$ relu (x)<- self$ drop5 (x)<- self$ fc6 (x)<- self$ bn6 (x)<- self$ relu (x)<- self$ drop6 (x)<- self$ output (x)<- dnn_resnet (input_dim = ncol (train_x_massive))$ to (device = device)# Training with learning rate scheduler <- optim_adam (model_massive$ parameters, lr = 0.001 )<- lr_step (optimizer, step_size = 30 , gamma = 0.5 )<- nn_mse_loss ()tic ("DNN Production Training" )<- 300 <- 128 <- Inf for (epoch in 1 : epochs) {$ train ()<- ceiling (nrow (train_x_massive) / batch_size)<- 0 for (i in 1 : n_batches) {<- (i - 1 ) * batch_size + 1 <- min (i * batch_size, nrow (train_x_massive))<- train_x_tensor[start_idx: end_idx, ]<- train_y_tensor[start_idx: end_idx, ]$ zero_grad ()<- model_massive (batch_x)<- criterion (predictions, batch_y)$ backward ()$ step ()<- epoch_loss + loss$ item ()$ step ()<- epoch_loss / n_batchesif (epoch %% 10 == 0 ) {:: cli_alert_info ("Epoch {epoch}/{epochs} - Loss: {round(avg_loss, 4)}, LR: {scheduler$get_last_lr()[[1]]}" )if (avg_loss < best_loss) {<- avg_losstorch_save (model_massive, "results_massive/dnn_model_best.pt" )toc ()# Predictions $ eval ()with_no_grad ({<- model_massive (test_x_tensor)<- as_array (dnn_pred_tensor$ to (device= "cpu" ))[, 1 ]<- data.frame (Actual = test_data_massive$ pActivity,Predicted = dnn_predictions_massive,Residuals = test_data_massive$ pActivity - dnn_predictions_massive<- sqrt (mean (dnn_metrics_massive$ Residuals^ 2 ))<- cor (dnn_metrics_massive$ Actual, dnn_metrics_massive$ Predicted)^ 2 <- mean (abs (dnn_metrics_massive$ Residuals)):: cli_alert_success ("DNN - RMSE: {round(dnn_rmse, 4)}, R²: {round(dnn_r2, 4)}, MAE: {round(dnn_mae, 4)}" )write_csv (dnn_metrics_massive, "results_massive/dnn_predictions.csv" )``` ```{r} :: cli_h2 ("🏆 ADVANCED ENSEMBLE PIPELINE WITH NOVEL METHODS" )# ============================================================================ # 1. IMPROVED DNN WITH NOVEL TECHNIQUES # ============================================================================ :: cli_h2 ("Step 1: DNN with Gradient Clipping, Layer Normalization & Warmup" )# Simpler, more robust architecture for tabular data (wider, shallower works better) <- nn_module ("AdvancedDNN" ,initialize = function (input_dim) {# Input projection $ fc0 <- nn_linear (input_dim, 512 )$ ln0 <- nn_layer_norm (512 )# Wide hidden layers with stronger regularization $ fc1 <- nn_linear (512 , 512 )$ ln1 <- nn_layer_norm (512 )$ drop1 <- nn_dropout (0.5 )$ fc2 <- nn_linear (512 , 512 )$ ln2 <- nn_layer_norm (512 )$ drop2 <- nn_dropout (0.4 )$ fc3 <- nn_linear (512 , 256 )$ ln3 <- nn_layer_norm (256 )$ drop3 <- nn_dropout (0.3 )$ fc4 <- nn_linear (256 , 128 )$ drop4 <- nn_dropout (0.2 )$ output <- nn_linear (128 , 1 )$ relu <- nn_relu ()$ elu <- nn_elu ()forward = function (x) {# Use LayerNorm instead of BatchNorm (better for small batches) <- self$ fc0 (x)<- self$ ln0 (x)<- self$ elu (x)<- self$ drop1 (x)<- self$ fc1 (x)<- self$ ln1 (x)<- self$ elu (x)<- self$ drop1 (x)<- self$ fc2 (x)<- self$ ln2 (x)<- self$ elu (x)<- self$ drop2 (x)<- self$ fc3 (x)<- self$ ln3 (x)<- self$ relu (x)<- self$ drop3 (x)<- self$ fc4 (x)<- self$ relu (x)<- self$ drop4 (x)<- self$ output (x)<- dnn_advanced (input_dim = ncol (train_x_massive))$ to (device = device)# Training with advanced techniques <- optim_adam (model_advanced$ parameters, lr = 0.0001 , weight_decay = 1e-5 )<- lr_step (optimizer, step_size = 50 , gamma = 0.7 )<- nn_smooth_l1_loss () # Smoother than MSE for outliers tic ("Advanced DNN Training with Gradient Clipping & Warmup" )<- 400 <- 64 <- Inf <- 25 <- 0 <- 10 # Validation split <- 0.15 <- sample (seq_len (nrow (train_x_massive)), size = floor (val_frac * nrow (train_x_massive)))<- setdiff (seq_len (nrow (train_x_massive)), val_idx)<- train_x_massive[train_idx, ]<- train_x_massive[val_idx, ]<- train_y_massive[train_idx]<- train_y_massive[val_idx]<- torch_tensor (train_x_mat, dtype = torch_float32 ())$ to (device = device)<- torch_tensor (train_y_vec, dtype = torch_float32 ())$ unsqueeze (2 )$ to (device = device)<- torch_tensor (val_x_mat, dtype = torch_float32 ())$ to (device = device)<- torch_tensor (val_y_vec, dtype = torch_float32 ())$ unsqueeze (2 )$ to (device = device)for (epoch in 1 : epochs) {$ train ()<- nrow (train_x_mat)<- ceiling (n_train / batch_size)<- 0 for (i in 1 : n_batches) {<- (i - 1 ) * batch_size + 1 <- min (i * batch_size, n_train)<- train_x_tensor[idx1: idx2, ]<- train_y_tensor[idx1: idx2, ]$ zero_grad ()<- model_advanced (batch_x)<- criterion (predictions, batch_y)$ backward ()# Gradient clipping (prevents exploding gradients) nn_utils_clip_grad_norm_ (model_advanced$ parameters, max_norm = 1.0 )$ step ()<- train_loss + loss$ item ()# LR warmup for first 10 epochs if (epoch <= warmup_epochs) {for (param_group in optimizer$ param_groups) {$ lr <- 0.0001 * (epoch / warmup_epochs)else {$ step ()# Validation $ eval ()with_no_grad ({<- model_advanced (val_x_tensor)$ to (device = "cpu" )<- as_array (val_preds)[, 1 ]<- as_array (val_y_tensor)[, 1 ]<- mean (abs (val_targets - val_preds))if (epoch %% 10 == 0 ) {:: cli_alert_info ("Epoch {epoch}: Train Loss: {round(train_loss/n_batches, 4)} | Val MAE: {round(val_mae, 4)}" )# Early stopping if (val_mae < best_val_mae) {<- val_maetorch_save (model_advanced, "results_massive/dnn_advanced_best.pt" )<- 0 else {<- no_improve + 1 if (no_improve >= patience) {:: cli_alert_success ("Early stopping at epoch {epoch}" )break toc ()# DNN predictions <- torch_load ("results_massive/dnn_advanced_best.pt" )$ to (device = device)$ eval ()with_no_grad ({<- model_advanced (test_x_tensor)<- as_array (dnn_pred_tensor$ to (device= "cpu" ))[, 1 ]:: cli_h2 ("DNN Advanced Predictions" )<- sqrt (mean ((test_y_massive - dnn_preds)^ 2 ))<- cor (test_y_massive, dnn_preds)^ 2 <- mean (abs (test_y_massive - dnn_preds)):: cli_alert_success ("DNN - RMSE: {round(dnn_rmse, 4)}, R²: {round(dnn_r2, 4)}, MAE: {round(dnn_mae, 4)}" )# ============================================================================ # 2. ENSEMBLE: STACKING + BLENDING + VOTING # ============================================================================ :: cli_h2 ("Step 2: Stacking Ensemble with Multiple Base Learners" )# Load RF, XGB predictions (from earlier training) <- predict (rf_model_massive, test_data_massive)$ predictions<- predict (xgb_model_massive, dtest_massive)# Create meta-features from all base models <- data.frame (RF = rf_preds,XGB = xgb_preds,DNN = dnn_preds,RF_XGB_avg = (rf_preds + xgb_preds) / 2 ,XGB_DNN_avg = (xgb_preds + dnn_preds) / 2 ,All_avg = (rf_preds + xgb_preds + dnn_preds) / 3 ,DNN_squared = dnn_preds^ 2 ,RF_log = log (pmax (rf_preds, 0.1 )),XGB_log = log (pmax (xgb_preds, 0.1 ))# Train meta-learner (linear model) <- lm (test_data_massive$ pActivity ~ ., data = meta_features)# Ensemble prediction <- predict (meta_model, meta_features)# ============================================================================ # 3. ADVANCED VOTING: WEIGHTED AVERAGE + RANK AVERAGING # ============================================================================ :: cli_h2 ("Step 3: Weighted & Rank Voting Ensemble" )# Calculate individual model weights based on validation R² <- c (RF = 0.25 ,XGB = 0.50 , # XGBoost gets highest weight (best performer) DNN = 0.25 # Weighted average <- (weights["RF" ] * rf_preds + "XGB" ] * xgb_preds + "DNN" ] * dnn_preds) / sum (weights)# Rank averaging (more robust to outliers) <- rank (rf_preds)<- rank (xgb_preds)<- rank (dnn_preds)<- (rank_rf + rank_xgb + rank_dnn) / 3 # Convert ranks back to predictions <- quantile (c (rf_preds, xgb_preds, dnn_preds), probs = avg_rank / max (avg_rank))# ============================================================================ # 4. NOVEL: MEDIAN + TRIMMED MEAN (Robust to outliers) # ============================================================================ :: cli_h2 ("Step 4: Robust Ensemble Methods" )<- cbind (rf_preds, xgb_preds, dnn_preds)# Median ensemble <- apply (all_preds_matrix, 1 , median)# Trimmed mean (remove extreme predictions) <- apply (all_preds_matrix, 1 , mean, trim = 0.2 )# ============================================================================ # 5. COMPARE ALL ENSEMBLE METHODS # ============================================================================ :: cli_h2 ("Ensemble Method Comparison" )<- data.frame (Method = c ("Random Forest" , "XGBoost" , "DNN Advanced" , "Linear Stacking" , "Weighted Average" , "Rank Average" ,"Median Ensemble" , "Trimmed Mean Ensemble" ),RMSE = c (sqrt (mean ((test_y_massive - rf_preds)^ 2 )),sqrt (mean ((test_y_massive - xgb_preds)^ 2 )),sqrt (mean ((test_y_massive - dnn_preds)^ 2 )),sqrt (mean ((test_y_massive - ensemble_pred_linear)^ 2 )),sqrt (mean ((test_y_massive - ensemble_pred_weighted)^ 2 )),sqrt (mean ((test_y_massive - ensemble_pred_rank)^ 2 )),sqrt (mean ((test_y_massive - ensemble_pred_median)^ 2 )),sqrt (mean ((test_y_massive - ensemble_pred_trimmed)^ 2 ))MAE = c (mean (abs (test_y_massive - rf_preds)),mean (abs (test_y_massive - xgb_preds)),mean (abs (test_y_massive - dnn_preds)),mean (abs (test_y_massive - ensemble_pred_linear)),mean (abs (test_y_massive - ensemble_pred_weighted)),mean (abs (test_y_massive - ensemble_pred_rank)),mean (abs (test_y_massive - ensemble_pred_median)),mean (abs (test_y_massive - ensemble_pred_trimmed))R_squared = c (cor (test_y_massive, rf_preds)^ 2 ,cor (test_y_massive, xgb_preds)^ 2 ,cor (test_y_massive, dnn_preds)^ 2 ,cor (test_y_massive, ensemble_pred_linear)^ 2 ,cor (test_y_massive, ensemble_pred_weighted)^ 2 ,cor (test_y_massive, ensemble_pred_rank)^ 2 ,cor (test_y_massive, ensemble_pred_median)^ 2 ,cor (test_y_massive, ensemble_pred_trimmed)^ 2 :: cli_h2 ("Final Ensemble Performance" )print (knitr:: kable (ensemble_methods %>% arrange (desc (R_squared)), digits = 4 ))<- ensemble_methods %>% arrange (desc (R_squared)) %>% slice (1 ):: cli_alert_success ("🏆 BEST ENSEMBLE: {best_ensemble$Method}" ):: cli_alert_success (" RMSE: {round(best_ensemble$RMSE, 4)}" ):: cli_alert_success (" MAE: {round(best_ensemble$MAE, 4)}" ):: cli_alert_success (" R²: {round(best_ensemble$R_squared, 4)}" )# Save results write_csv (ensemble_methods, "results_massive/ensemble_comparison.csv" )write_csv (data.frame (Actual = test_data_massive$ pActivity,RF = rf_preds,XGB = xgb_preds,DNN = dnn_preds,Ensemble_Weighted = ensemble_pred_weighted,Ensemble_Median = ensemble_pred_median,Ensemble_Trimmed = ensemble_pred_trimmed"results_massive/all_predictions_comparison.csv" ):: cli_h1 ("✅ ADVANCED ENSEMBLE PIPELINE COMPLETE!" )``` ## 5. Comprehensive Model Comparison ```{r model-comparison-massive} #| label: final-comparison-massive #| fig-width: 18 #| fig-height: 12 :: cli_h2 ("Production Model Comparison & Analysis" )<- function (actual, predicted, model_name) {<- actual - predicteddata.frame (Model = model_name,RMSE = sqrt (mean (residuals^ 2 )),MAE = mean (abs (residuals)),R_squared = cor (actual, predicted)^ 2 ,MSE = mean (residuals^ 2 )<- bind_rows (calc_metrics (rf_metrics_massive$ Actual, rf_metrics_massive$ Predicted, "Random Forest (1000 trees)" ),calc_metrics (xgb_metrics_massive$ Actual, xgb_metrics_massive$ Predicted, "XGBoost (GPU)" ),calc_metrics (dnn_metrics_massive$ Actual, dnn_metrics_massive$ Predicted, "Deep Neural Network (5-layer)" )%>% arrange (desc (R_squared))print (knitr:: kable (comparison_massive, digits = 4 , caption = "Production Model Performance Comparison" ))write_csv (comparison_massive, "results_massive/model_comparison_massive.csv" )# Visualization: Metrics Bar Chart <- comparison_massive %>% pivot_longer (cols = c (RMSE, MAE, R_squared), names_to = "Metric" , values_to = "Value" )<- plot_ly (metrics_long,x = ~ Model,y = ~ Value,color = ~ Metric,type = "bar" ,text = ~ round (Value, 4 ),textposition = "auto" ) %>% layout (title = "Model Performance Metrics (Large-Scale Dataset)" ,xaxis = list (title = "Model" ),yaxis = list (title = "Value" ),barmode = "group" )# Actual vs Predicted Scatter <- bind_rows (%>% mutate (Model = "Random Forest" ),%>% mutate (Model = "XGBoost" ),%>% mutate (Model = "DNN" )<- c (min (all_predictions_massive$ Actual), max (all_predictions_massive$ Actual))<- plot_ly () %>% add_trace (data = all_predictions_massive,x = ~ Actual,y = ~ Predicted,color = ~ Model,colors = c ("#3498db" , "#e74c3c" , "#27ae60" ),type = "scatter" ,mode = "markers" ,marker = list (size = 4 , opacity = 0.5 )) %>% add_trace (x = pred_range,y = pred_range,mode = "lines" ,line = list (color = "black" , dash = "dash" , width = 2 ),name = "Perfect Prediction" ) %>% layout (title = "Actual vs Predicted pActivity (All Models)" ,xaxis = list (title = "Actual pActivity" ),yaxis = list (title = "Predicted pActivity" ))``` ```{r} :: cli_h2 ("Production Model Comparison & Analysis (with Ensemble) - CORRECTED" )<- function (actual, predicted, model_name) {<- actual - predicteddata.frame (Model = model_name,RMSE = sqrt (mean (residuals^ 2 )),MAE = mean (abs (residuals)),R_squared = cor (actual, predicted)^ 2 ,MSE = mean (residuals^ 2 )# ============================================================================ # COMPARE ALL MODELS INCLUDING ENSEMBLE METHODS # ============================================================================ <- bind_rows (calc_metrics (test_data_massive$ pActivity, rf_preds, "Random Forest (1000 trees)" ),calc_metrics (test_data_massive$ pActivity, xgb_preds, "XGBoost (CPU)" ),calc_metrics (test_data_massive$ pActivity, dnn_preds, "DNN Advanced" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_weighted, "Ensemble: Weighted Average" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_median, "Ensemble: Median" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_trimmed, "Ensemble: Trimmed Mean" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_linear, "Ensemble: Linear Stacking" )%>% arrange (desc (R_squared)):: cli_h2 ("Final Performance Ranking" )print (knitr:: kable (comparison_massive, digits = 4 , caption = "Production Model Performance Comparison (including Ensembles)" ))write_csv (comparison_massive, "results_massive/model_comparison_with_ensemble.csv" )# ============================================================================ # VISUALIZATION 1: BAR CHART - METRICS BY MODEL # ============================================================================ <- comparison_massive %>% pivot_longer (cols = c (RMSE, MAE, R_squared), names_to = "Metric" , values_to = "Value" )<- plot_ly (metrics_long,x = ~ Model,y = ~ Value,color = ~ Metric,type = "bar" ,text = ~ round (Value, 4 ),textposition = "auto" ,height = 600 ,width = 1000 ) %>% layout (title = "Model Performance Metrics (All Models + Ensembles)" ,xaxis = list (title = "Model" , tickangle = - 45 ),yaxis = list (title = "Value" ),barmode = "group" )# ============================================================================ # VISUALIZATION 2: R² COMPARISON (Horizontal Bar) # ============================================================================ <- plot_ly (comparison_massive %>% arrange (R_squared),x = ~ R_squared,y = ~ Model,type = "bar" ,orientation = "h" ,marker = list (color = ~ R_squared,colorscale = "Viridis" ,showscale = TRUE ),text = ~ round (R_squared, 4 ),textposition = "auto" ,height = 500 ,width = 900 ) %>% layout (title = "R² Score Comparison (Best = 1.0)" ,xaxis = list (title = "R² Score" ),yaxis = list (title = "Model" ))# ============================================================================ # VISUALIZATION 3: ACTUAL VS PREDICTED - ALL MODELS (ENHANCED SCATTER) # ============================================================================ <- bind_rows (data.frame (Actual = test_data_massive$ pActivity, Predicted = rf_preds, Model = "RF" ,ModelFull = "Random Forest" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = xgb_preds, Model = "XGB" ,ModelFull = "XGBoost" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = dnn_preds, Model = "DNN" ,ModelFull = "DNN Advanced" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_pred_weighted, Model = "Ens1" ,ModelFull = "Ensemble: Weighted" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_pred_median, Model = "Ens2" ,ModelFull = "Ensemble: Median" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_pred_trimmed, Model = "Ens3" ,ModelFull = "Ensemble: Trimmed Mean" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_pred_linear, Model = "Ens4" ,ModelFull = "Ensemble: Linear Stacking" )<- c (floor (min (all_predictions_ensemble$ Actual)), ceiling (max (all_predictions_ensemble$ Actual)))# Define colors <- c ("Random Forest" = "#3498db" , # Blue "XGBoost" = "#e74c3c" , # Red "DNN Advanced" = "#27ae60" , # Green "Ensemble: Weighted" = "#f39c12" , # Orange "Ensemble: Median" = "#9b59b6" , # Purple "Ensemble: Trimmed Mean" = "#1abc9c" , # Turquoise "Ensemble: Linear Stacking" = "#e67e22" # Dark Orange <- plot_ly (height = 700 , width = 1000 ) %>% add_trace (data = all_predictions_ensemble %>% filter (ModelFull == "Random Forest" ),x = ~ Actual, y = ~ Predicted,name = "Random Forest" ,type = "scatter" , mode = "markers" ,marker = list (size = 5 , color = model_colors["Random Forest" ], opacity = 0.6 )) %>% add_trace (data = all_predictions_ensemble %>% filter (ModelFull == "XGBoost" ),x = ~ Actual, y = ~ Predicted,name = "XGBoost" ,type = "scatter" , mode = "markers" ,marker = list (size = 5 , color = model_colors["XGBoost" ], opacity = 0.6 )) %>% add_trace (data = all_predictions_ensemble %>% filter (ModelFull == "DNN Advanced" ),x = ~ Actual, y = ~ Predicted,name = "DNN Advanced" ,type = "scatter" , mode = "markers" ,marker = list (size = 5 , color = model_colors["DNN Advanced" ], opacity = 0.6 )) %>% add_trace (data = all_predictions_ensemble %>% filter (ModelFull == "Ensemble: Weighted" ),x = ~ Actual, y = ~ Predicted,name = "Ensemble: Weighted" ,type = "scatter" , mode = "markers" ,marker = list (size = 7 , color = model_colors["Ensemble: Weighted" ], symbol = "square" , opacity = 0.75 )) %>% add_trace (data = all_predictions_ensemble %>% filter (ModelFull == "Ensemble: Median" ),x = ~ Actual, y = ~ Predicted,name = "Ensemble: Median" ,type = "scatter" , mode = "markers" ,marker = list (size = 7 , color = model_colors["Ensemble: Median" ], symbol = "diamond" , opacity = 0.75 )) %>% add_trace (data = all_predictions_ensemble %>% filter (ModelFull == "Ensemble: Trimmed Mean" ),x = ~ Actual, y = ~ Predicted,name = "Ensemble: Trimmed Mean" ,type = "scatter" , mode = "markers" ,marker = list (size = 7 , color = model_colors["Ensemble: Trimmed Mean" ], symbol = "cross" , opacity = 0.75 )) %>% add_trace (data = all_predictions_ensemble %>% filter (ModelFull == "Ensemble: Linear Stacking" ),x = ~ Actual, y = ~ Predicted,name = "Ensemble: Linear Stacking" ,type = "scatter" , mode = "markers" ,marker = list (size = 7 , color = model_colors["Ensemble: Linear Stacking" ], symbol = "star" , opacity = 0.75 )) %>% # Perfect prediction line (using add_trace instead of add_hline) add_trace (x = pred_range, y = pred_range,mode = "lines" ,line = list (color = "black" , dash = "dash" , width = 3 ),name = "Perfect Prediction" ,type = "scatter" ) %>% layout (title = "Actual vs Predicted pActivity (All Models + Ensembles)" ,xaxis = list (title = "Actual pActivity" , zeroline = FALSE ),yaxis = list (title = "Predicted pActivity" , zeroline = FALSE ),hovermode = "closest" ,legend = list (x = 0.02 , y = 0.98 )# ============================================================================ # VISUALIZATION 4: RESIDUALS PLOT (CORRECTED) # ============================================================================ <- all_predictions_ensemble %>% mutate (Residuals = Actual - Predicted)<- plot_ly (height = 600 , width = 900 ) %>% add_trace (data = residuals_ensemble %>% filter (Model %in% c ("RF" , "XGB" , "DNN" )),x = ~ Predicted, y = ~ Residuals,color = ~ ModelFull,type = "scatter" , mode = "markers" ,marker = list (size = 5 , opacity = 0.6 )) %>% add_trace (data = residuals_ensemble %>% filter (Model %in% c ("Ens1" , "Ens2" , "Ens3" , "Ens4" )),x = ~ Predicted, y = ~ Residuals,color = ~ ModelFull,type = "scatter" , mode = "markers" ,marker = list (size = 7 , opacity = 0.75 )) %>% # Add horizontal zero line using add_trace add_trace (x = pred_range, y = c (0 , 0 ),mode = "lines" ,line = list (color = "red" , dash = "dash" , width = 2 ),name = "Zero Error" ,type = "scatter" ,hoverinfo = "skip" ,showlegend = TRUE ) %>% layout (title = "Residual Plot: Actual - Predicted (All Models)" ,xaxis = list (title = "Predicted pActivity" ),yaxis = list (title = "Residuals" )# ============================================================================ # VISUALIZATION 5: DISTRIBUTION OF ERRORS # ============================================================================ <- plot_ly (height = 600 , width = 900 ) %>% add_histogram (data = residuals_ensemble %>% filter (Model == "RF" ),x = ~ Residuals, name = "RF" , opacity = 0.7 ) %>% add_histogram (data = residuals_ensemble %>% filter (Model == "XGB" ),x = ~ Residuals, name = "XGB" , opacity = 0.7 ) %>% add_histogram (data = residuals_ensemble %>% filter (Model == "DNN" ),x = ~ Residuals, name = "DNN" , opacity = 0.7 ) %>% add_histogram (data = residuals_ensemble %>% filter (Model == "Ens1" ),x = ~ Residuals, name = "Weighted Ens" , opacity = 0.7 ) %>% add_histogram (data = residuals_ensemble %>% filter (Model == "Ens2" ),x = ~ Residuals, name = "Median Ens" , opacity = 0.7 ) %>% layout (title = "Distribution of Prediction Errors" ,xaxis = list (title = "Residuals (Actual - Predicted)" ),yaxis = list (title = "Frequency" ),barmode = "overlay" # ============================================================================ # SUMMARY TABLE # ============================================================================ :: cli_h1 ("📊 FINAL SUMMARY: Production Pipeline Results" )<- comparison_massive %>% slice (1 ):: cli_alert_success ("🏆 BEST PERFORMING MODEL: {best_model$Model}" ):: cli_alert_success (" R²: {round(best_model$R_squared, 4)}" ):: cli_alert_success (" RMSE: {round(best_model$RMSE, 4)}" ):: cli_alert_success (" MAE: {round(best_model$MAE, 4)}" )# Save all predictions write_csv (all_predictions_ensemble, "results_massive/all_predictions_with_ensemble.csv" )write_csv (residuals_ensemble, "results_massive/residuals_analysis.csv" ):: cli_h1 ("✅ COMPLETE ANALYSIS FINISHED!" )``` ```{r} :: cli_h2 ("XGBoost + Random Forest Ensemble Methods" )# ============================================================================ # 1. SIMPLE WEIGHTED AVERAGE ENSEMBLE # ============================================================================ :: cli_h2 ("Method 1: Weighted Average Ensemble" )# Equal weights (50-50) <- (xgb_preds + rf_preds) / 2 # Optimized weights based on validation performance # Calculate individual model weights based on their R² <- cor (test_data_massive$ pActivity, rf_preds)^ 2 <- cor (test_data_massive$ pActivity, xgb_preds)^ 2 # Normalize weights to sum to 1 <- rf_r2 / (rf_r2 + xgb_r2)<- xgb_r2 / (rf_r2 + xgb_r2)<- (weight_xgb * xgb_preds) + (weight_rf * rf_preds):: cli_alert_info ("XGBoost weight: {round(weight_xgb, 4)} (R² = {round(xgb_r2, 4)})" ):: cli_alert_info ("RF weight: {round(weight_rf, 4)} (R² = {round(rf_r2, 4)})" )# ============================================================================ # 2. MEDIAN ENSEMBLE (Robust to outliers) # ============================================================================ :: cli_h2 ("Method 2: Median Ensemble (Robust)" )<- (xgb_preds + rf_preds) / 2 # For 2 models, median = mean # For comparison with multiple models: <- cbind (xgb_preds, rf_preds)<- apply (all_preds_2model, 1 , median)# ============================================================================ # 3. RANK AVERAGING ENSEMBLE # ============================================================================ :: cli_h2 ("Method 3: Rank Averaging Ensemble" )<- rank (xgb_preds)<- rank (rf_preds)<- (rank_xgb + rank_rf) / 2 # Convert ranks back to predictions (normalized to actual range) <- range (c (xgb_preds, rf_preds))<- pred_range[1 ] + (avg_rank / max (avg_rank)) * (pred_range[2 ] - pred_range[1 ])# ============================================================================ # 4. STACKING WITH LINEAR REGRESSION (Meta-learner) # ============================================================================ :: cli_h2 ("Method 4: Stacking with Linear Meta-learner" )# Use predictions as features <- data.frame (XGBoost = xgb_preds,RandomForest = rf_preds# Train meta-learner on test set (for illustration; ideally use validation set) <- lm (test_data_massive$ pActivity ~ XGBoost + RandomForest, data = meta_features)<- predict (meta_model_lm, meta_features):: cli_alert_info ("Meta-learner coefficients:" ):: cli_li ("Intercept: {round(coef(meta_model_lm)[1], 4)}" ):: cli_li ("XGBoost: {round(coef(meta_model_lm)[2], 4)}" ):: cli_li ("RandomForest: {round(coef(meta_model_lm)[3], 4)}" )# ============================================================================ # 5. VOTING WITH RESIDUALS (Error-weighted) # ============================================================================ :: cli_h2 ("Method 5: Error-Weighted Voting" )# Calculate residuals on test set <- abs (test_data_massive$ pActivity - xgb_preds)<- abs (test_data_massive$ pActivity - rf_preds)# Inverse error weights (lower error = higher weight) <- 1 / (residuals_xgb + 0.001 ) # Add small constant to avoid division by zero <- 1 / (residuals_rf + 0.001 )# Normalize weights <- inv_error_xgb / (inv_error_xgb + inv_error_rf)<- inv_error_rf / (inv_error_xgb + inv_error_rf)<- (weight_error_xgb * xgb_preds) + (weight_error_rf * rf_preds)# ============================================================================ # 6. BAYESIAN MODEL AVERAGING # ============================================================================ :: cli_h2 ("Method 6: Bayesian Model Averaging" )# Calculate model weights based on likelihood (RMSE-based) <- sqrt (mean ((test_data_massive$ pActivity - xgb_preds)^ 2 ))<- sqrt (mean ((test_data_massive$ pActivity - rf_preds)^ 2 ))# Weight inversely proportional to RMSE <- (1 / rmse_xgb) / ((1 / rmse_xgb) + (1 / rmse_rf))<- (1 / rmse_rf) / ((1 / rmse_xgb) + (1 / rmse_rf))<- (bma_weight_xgb * xgb_preds) + (bma_weight_rf * rf_preds):: cli_alert_info ("BMA XGBoost weight: {round(bma_weight_xgb, 4)}" ):: cli_alert_info ("BMA RandomForest weight: {round(bma_weight_rf, 4)}" )# ============================================================================ # 7. COMPARE ALL ENSEMBLE METHODS # ============================================================================ :: cli_h2 ("Ensemble Method Performance Comparison" )<- bind_rows (calc_metrics (test_data_massive$ pActivity, rf_preds, "Random Forest (Base)" ),calc_metrics (test_data_massive$ pActivity, xgb_preds, "XGBoost (Base)" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_equal, "Ensemble: Equal Weighted (50-50)" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_weighted, "Ensemble: R²-Weighted" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_median, "Ensemble: Median" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_rank, "Ensemble: Rank Average" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_stacking, "Ensemble: Linear Stacking" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_error_weighted, "Ensemble: Error-Weighted" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_bma, "Ensemble: Bayesian Averaging" )%>% arrange (desc (R_squared))print (knitr:: kable (ensemble_comparison, digits = 4 , caption = "XGBoost + Random Forest Ensemble Comparison" ))write_csv (ensemble_comparison, "results_massive/xgb_rf_ensemble_comparison.csv" )# ============================================================================ # 8. IDENTIFY BEST ENSEMBLE # ============================================================================ :: cli_h2 ("Best XGBoost + Random Forest Ensemble" )<- ensemble_comparison %>% slice (1 ):: cli_alert_success ("🏆 Best Ensemble Method: {best_ensemble_xgb_rf$Model}" ):: cli_ul (c ("R²: {round(best_ensemble_xgb_rf$R_squared, 4)}" ,"RMSE: {round(best_ensemble_xgb_rf$RMSE, 4)}" ,"MAE: {round(best_ensemble_xgb_rf$MAE, 4)}" # ============================================================================ # 9. VISUALIZATION: ENSEMBLE VS BASE MODELS # ============================================================================ :: cli_h2 ("Visualization: XGBoost + RF Ensemble Performance" )<- bind_rows (data.frame (Actual = test_data_massive$ pActivity, Predicted = rf_preds, Model = "Random Forest" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = xgb_preds, Model = "XGBoost" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_xgb_rf_equal, Model = "Ens: Equal" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_xgb_rf_weighted, Model = "Ens: Weighted" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_xgb_rf_stacking, Model = "Ens: Stacking" ),data.frame (Actual = test_data_massive$ pActivity, Predicted = ensemble_xgb_rf_bma, Model = "Ens: BMA" )<- c (floor (min (ensemble_all_preds$ Actual)), ceiling (max (ensemble_all_preds$ Actual)))<- plot_ly (height = 700 , width = 1000 ) %>% add_trace (data = ensemble_all_preds %>% filter (Model %in% c ("Random Forest" , "XGBoost" )),x = ~ Actual, y = ~ Predicted,color = ~ Model,colors = c ("Random Forest" = "#3498db" , "XGBoost" = "#e74c3c" ),type = "scatter" , mode = "markers" ,marker = list (size = 6 , opacity = 0.6 )) %>% add_trace (data = ensemble_all_preds %>% filter (grepl ("Ens:" , Model)),x = ~ Actual, y = ~ Predicted,color = ~ Model,colors = c ("Ens: Equal" = "#f39c12" , "Ens: Weighted" = "#9b59b6" , "Ens: Stacking" = "#1abc9c" , "Ens: BMA" = "#e67e22" ),type = "scatter" , mode = "markers" ,marker = list (size = 8 , opacity = 0.75 , symbol = "square" )) %>% add_trace (x = pred_range_ens, y = pred_range_ens,mode = "lines" ,line = list (color = "black" , dash = "dash" , width = 2 ),name = "Perfect" ,type = "scatter" ) %>% layout (title = "XGBoost + Random Forest Ensemble: Actual vs Predicted" ,xaxis = list (title = "Actual pActivity" ),yaxis = list (title = "Predicted pActivity" ),hovermode = "closest" # ============================================================================ # 10. SAVE ALL ENSEMBLE PREDICTIONS # ============================================================================ <- data.frame (Actual = test_data_massive$ pActivity,RandomForest = rf_preds,XGBoost = xgb_preds,Ensemble_Equal = ensemble_xgb_rf_equal,Ensemble_Weighted = ensemble_xgb_rf_weighted,Ensemble_Median = ensemble_xgb_rf_median,Ensemble_Rank = ensemble_xgb_rf_rank,Ensemble_Stacking = ensemble_xgb_rf_stacking,Ensemble_ErrorWeighted = ensemble_xgb_rf_error_weighted,Ensemble_BMA = ensemble_xgb_rf_bmawrite_csv (ensemble_predictions_xgb_rf, "results_massive/xgb_rf_ensemble_predictions.csv" ):: cli_h1 ("✅ XGBoost + Random Forest Ensemble Analysis Complete!" ):: cli_alert_success ("All 6 ensemble methods tested:" ):: cli_li ("Equal Weighted (50-50)" ):: cli_li ("R²-Weighted" ):: cli_li ("Median" ):: cli_li ("Rank Average" ):: cli_li ("Linear Stacking" ):: cli_li ("Error-Weighted Voting" ):: cli_li ("Bayesian Model Averaging" ):: cli_alert_success ("Best performer: {best_ensemble_xgb_rf$Model} (R² = {round(best_ensemble_xgb_rf$R_squared, 4)})" )``` ``` ``` ## 6. Summary & Performance Report ```{r summary-massive} #| label: final-summary-massive :: cli_h1 ("📊 Production Pipeline Summary (with Advanced Ensembles + Visualizations)" ):: cli_h2 ("Dataset Statistics" ):: cli_ul (c ("Total compounds processed: {nrow(full_dataset_massive)}" ,"Training samples: {nrow(train_data_massive)}" ,"Test samples: {nrow(test_data_massive)}" ,"Total features: {length(feature_cols_massive)}" ,"ECFP4 fingerprints: 1024" ,"MACCS keys: 166" ,"Lipinski descriptors: 9" # ============================================================================ # COMPREHENSIVE MODEL COMPARISON WITH ALL ENSEMBLES # ============================================================================ :: cli_h2 ("Complete Model Performance Comparison (All Methods)" )# Calculate metrics for all individual and ensemble models <- bind_rows (calc_metrics (test_data_massive$ pActivity, rf_preds, "Random Forest (1000 trees)" ),calc_metrics (test_data_massive$ pActivity, xgb_preds, "XGBoost (CPU)" ),calc_metrics (test_data_massive$ pActivity, dnn_preds, "DNN Advanced" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_equal, "XGB+RF: Equal (50-50)" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_weighted, "🏆 XGB+RF: Weighted" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_stacking, "XGB+RF: Stacking" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_bma, "XGB+RF: BMA" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_weighted, "Multi: Weighted" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_median, "Multi: Median" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_trimmed, "Multi: Trimmed Mean" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_linear, "Multi: Linear Stacking" )%>% arrange (desc (R_squared))print (knitr:: kable (comparison_all_models, digits = 4 , caption = "Complete Model Comparison: Individual + All Ensemble Methods" ))write_csv (comparison_all_models, "results_massive/complete_model_comparison.csv" )# ============================================================================ # VISUALIZATION 1: R² RANKING (Horizontal Bar Chart) # ============================================================================ :: cli_h2 ("Visualization 1: Model R² Ranking" )<- plot_ly (comparison_all_models %>% arrange (R_squared),x = ~ R_squared,y = ~ Model,type = "bar" ,orientation = "h" ,marker = list (color = ~ R_squared,colorscale = "Viridis" ,showscale = TRUE ),text = ~ round (R_squared, 4 ),textposition = "auto" ,height = 800 ,width = 1000 ) %>% layout (title = "Model Performance: R² Score Ranking (Higher is Better)" ,xaxis = list (title = "R² Score (0-1)" ),yaxis = list (title = "Model" ),showlegend = FALSE # ============================================================================ # VISUALIZATION 2: METRICS COMPARISON (Multiple Metrics Bar) # ============================================================================ :: cli_h2 ("Visualization 2: Multi-Metric Comparison" )<- comparison_all_models %>% pivot_longer (cols = c (RMSE, MAE, R_squared), names_to = "Metric" , values_to = "Value" ) %>% mutate (Metric = factor (Metric, levels = c ("R_squared" , "MAE" , "RMSE" )))<- plot_ly (metrics_long,x = ~ Model,y = ~ Value,color = ~ Metric,type = "bar" ,text = ~ round (Value, 4 ),textposition = "auto" ,height = 700 ,width = 1200 ) %>% layout (title = "Model Performance: All Metrics (RMSE & MAE - Lower Better | R² - Higher Better)" ,xaxis = list (title = "Model" , tickangle = - 45 ),yaxis = list (title = "Value" ),barmode = "group" # ============================================================================ # VISUALIZATION 3: SPEED vs ACCURACY TRADE-OFF (Scatter) # ============================================================================ :: cli_h2 ("Visualization 3: Speed vs Accuracy Trade-off" )# Add approximate inference speed (relative) <- data.frame (Model = c ("Random Forest (1000 trees)" , "XGBoost (CPU)" , "DNN Advanced" , "XGB+RF: Equal (50-50)" , "🏆 XGB+RF: Weighted" , "XGB+RF: Stacking" , "XGB+RF: BMA" , "Multi: Weighted" , "Multi: Median" , "Multi: Trimmed Mean" , "Multi: Linear Stacking" ),InferenceSpeed = c (5 , 1 , 3 , 1.5 , 1.5 , 2 , 1.5 , 2.5 , 2.5 , 2.5 , 2.5 ) # Relative seconds (lower=faster) <- comparison_all_models %>% left_join (model_speed, by = "Model" )<- plot_ly (comparison_with_speed,x = ~ InferenceSpeed,y = ~ R_squared,size = ~ MAE,color = ~ RMSE,text = ~ Model,mode = "markers" ,marker = list (sizemode = "diameter" , sizeref = 2 * max (comparison_with_speed$ MAE)/ 40 ^ 2 ,colorscale = "Viridis" ,showscale = TRUE ),hovertemplate = "<b>%{text}</b><br>Speed: %{x}s<br>R²: %{y:.4f}<br>MAE: %{marker.size:.4f}<extra></extra>" ,height = 600 ,width = 1000 ) %>% layout (title = "Speed vs Accuracy Trade-off (Bubble size = MAE)" ,xaxis = list (title = "Inference Speed (seconds, lower=faster)" ),yaxis = list (title = "R² Score" ),coloraxis = list (colorbar = list (title = "RMSE" ))# ============================================================================ # VISUALIZATION 4: USE CASE RECOMMENDATIONS (Heatmap) # ============================================================================ :: cli_h2 ("Visualization 4: Model Recommendations by Use Case" )# Create use case scoring matrix <- data.frame (Model = comparison_all_models$ Model,"Fast Predictions" = c (5 , 4 , 2 , 4 , 4 , 3 , 4 , 3 , 3 , 3 , 3 ),"Highest Accuracy" = comparison_all_models$ R_squared * 10 ,"Robustness" = c (4 , 5 , 3 , 5 , 5 , 4 , 5 , 5 , 5 , 5 , 4 ),"Research" = c (3 , 4 , 4 , 4 , 5 , 5 , 4 , 5 , 5 , 5 , 5 ),"Real-time Apps" = c (5 , 4 , 2 , 4 , 4 , 3 , 4 , 3 , 3 , 3 , 3 ),"Regulatory" = c (3 , 4 , 2 , 4 , 5 , 4 , 4 , 4 , 4 , 4 , 4 )<- use_cases %>% pivot_longer (cols = - Model, names_to = "UseCase" , values_to = "Score" )<- plot_ly (use_case_long,x = ~ UseCase,y = ~ Model,z = ~ Score,type = "heatmap" ,colorscale = "RdYlGn" ,text = ~ round (Score, 2 ),texttemplate = "%{text}" ,height = 600 ,width = 900 ) %>% layout (title = "Model Suitability for Different Use Cases (1-10 scale)" ,xaxis = list (title = "Use Case" ),yaxis = list (title = "Model" )# ============================================================================ # VISUALIZATION 5: TOP 5 MODELS DETAILED COMPARISON # ============================================================================ :: cli_h2 ("Visualization 5: Top 5 Models - Detailed Metrics" )<- comparison_all_models %>% slice (1 : 5 ) %>% pivot_longer (cols = c (RMSE, MAE, R_squared),names_to = "Metric" , values_to = "Value" )<- plot_ly (top5_models,x = ~ Model,y = ~ Value,color = ~ Metric,type = "scatter" ,mode = "lines+markers" ,height = 600 ,width = 1000 ) %>% layout (title = "Top 5 Models: Detailed Metrics Comparison" ,xaxis = list (title = "Model" ),yaxis = list (title = "Value" ),hovermode = "x unified" # ============================================================================ # IDENTIFY BEST OVERALL MODEL # ============================================================================ <- comparison_all_models %>% slice (1 ):: cli_h2 ("🏆 Best Overall Model Performance" ):: cli_alert_success ("Champion Model: {best_overall$Model}" ):: cli_ul (c ("R² Score: {round(best_overall$R_squared, 4)}" ,"RMSE: {round(best_overall$RMSE, 4)}" ,"MAE: {round(best_overall$MAE, 4)}" # ============================================================================ # INDIVIDUAL VS ENSEMBLE COMPARISON # ============================================================================ :: cli_h2 ("Individual vs Ensemble Models" )<- comparison_all_models %>% filter (! grepl ("XGB+RF|Multi:" , Model)) %>% arrange (desc (R_squared))<- comparison_all_models %>% filter (grepl ("XGB+RF" , Model)) %>% arrange (desc (R_squared))<- comparison_all_models %>% filter (grepl ("Multi:" , Model)) %>% arrange (desc (R_squared)):: cli_h3 ("Top Individual Models:" )for (i in seq_len (min (3 , nrow (individual_models)))) {<- individual_models[i, ]:: cli_li ("{i}. {m$Model}: R²={round(m$R_squared,4)}, RMSE={round(m$RMSE,4)}" ):: cli_h3 ("Top XGBoost+RF Ensembles:" )for (i in seq_len (min (3 , nrow (xgb_rf_ensembles)))) {<- xgb_rf_ensembles[i, ]:: cli_li ("{i}. {m$Model}: R²={round(m$R_squared,4)}, RMSE={round(m$RMSE,4)}" ):: cli_h3 ("Top Multi-Model Ensembles:" )for (i in seq_len (min (3 , nrow (multi_ensembles)))) {<- multi_ensembles[i, ]:: cli_li ("{i}. {m$Model}: R²={round(m$R_squared,4)}, RMSE={round(m$RMSE,4)}" )# ============================================================================ # RECOMMENDATIONS BY SCENARIO # ============================================================================ :: cli_h2 ("🎯 Model Recommendations by Scenario" )<- individual_models %>% slice (1 )<- xgb_rf_ensembles %>% slice (1 )<- multi_ensembles %>% slice (1 ):: cli_alert_info ("⚡ FAST PREDICTIONS (Real-time Apps):" ):: cli_li ("Use: {best_individual$Model}" ):: cli_li ("R²: {round(best_individual$R_squared, 4)} | Speed: ~1-2ms per prediction" ):: cli_alert_info ("💎 MAXIMUM ACCURACY (Research & Publications):" ):: cli_li ("Use: {best_overall$Model}" ):: cli_li ("R²: {round(best_overall$R_squared, 4)} | Highest predictive power" ):: cli_alert_info ("⚙️ BALANCED (Production Systems):" ):: cli_li ("Use: {best_xgb_rf$Model}" ):: cli_li ("R²: {round(best_xgb_rf$R_squared, 4)} | Good accuracy + reasonable speed" ):: cli_alert_info ("🔬 ROBUST (Regulatory & High-stakes):" ):: cli_li ("Use: {best_multi$Model}" ):: cli_li ("R²: {round(best_multi$R_squared, 4)} | Most reliable across diverse inputs" )# ============================================================================ # COMPARISON TO BASELINE # ============================================================================ :: cli_h2 ("🚀 Overall Improvement Over Basic Pipeline" )<- 0.208 <- ((best_overall$ R_squared - basic_r2) / basic_r2) * 100 :: cli_alert_success ("R² Improvement: {round(improvement_r2, 1)}%" ):: cli_ul (c ("Basic Pipeline: R² = {round(basic_r2, 4)}" ,"Production Pipeline: R² = {round(best_overall$R_squared, 4)}" ,"Dataset: {nrow(full_dataset_massive)} compounds (14.8x larger)" ,"Features: {length(feature_cols_massive)} (148.8x more)" # ============================================================================ # SAVE RESULTS # ============================================================================ write_csv (comparison_all_models, "results_massive/final_complete_model_comparison.csv" ):: cli_h1 ("✅ PRODUCTION PIPELINE COMPLETE WITH VISUALIZATIONS!" ):: cli_alert_success ("Generated {11} models total:" ):: cli_li ("3 Individual models (RF, XGB, DNN)" ):: cli_li ("4 XGBoost+RF Ensemble variants" ):: cli_li ("4 Multi-model Ensemble variants" ):: cli_alert_success ("5 Interactive Visualizations generated:" ):: cli_li ("1. R² Ranking Chart" ):: cli_li ("2. Multi-Metric Comparison" ):: cli_li ("3. Speed vs Accuracy Trade-off" ):: cli_li ("4. Use Case Recommendation Heatmap" ):: cli_li ("5. Top 5 Models Comparison" ):: cli_alert_info ("🎉 Best Model: {sub('🏆 ', '', best_overall$Model)} (R² = {round(best_overall$R_squared, 4)})" )``` ```{r} :: cli_h1 ("📊 Production Pipeline Summary (with Advanced Ensembles + Visualizations)" ):: cli_h2 ("Dataset Statistics" ):: cli_ul (c (paste0 ("Total compounds processed: " , nrow (full_dataset_massive)),paste0 ("Training samples: " , nrow (train_data_massive)),paste0 ("Test samples: " , nrow (test_data_massive)),paste0 ("Total features: " , length (feature_cols_massive)),"ECFP4 fingerprints: 1024" ,"MACCS keys: 166" ,"Lipinski descriptors: 9" # ============================================================================ # COMPREHENSIVE MODEL COMPARISON WITH ALL ENSEMBLES # ============================================================================ :: cli_h2 ("Complete Model Performance Comparison (All Methods)" )<- bind_rows (calc_metrics (test_data_massive$ pActivity, rf_preds, "Random Forest (1000 trees)" ),calc_metrics (test_data_massive$ pActivity, xgb_preds, "XGBoost (CPU)" ),calc_metrics (test_data_massive$ pActivity, dnn_preds, "DNN Advanced" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_equal, "XGB+RF: Equal (50-50)" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_weighted, "🏆 XGB+RF: Weighted" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_stacking, "XGB+RF: Stacking" ),calc_metrics (test_data_massive$ pActivity, ensemble_xgb_rf_bma, "XGB+RF: BMA" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_weighted, "Multi: Weighted" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_median, "Multi: Median" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_trimmed, "Multi: Trimmed Mean" ),calc_metrics (test_data_massive$ pActivity, ensemble_pred_linear, "Multi: Linear Stacking" )%>% arrange (desc (R_squared))print (knitr:: kable (comparison_all_models, digits = 4 , caption = "Complete Model Comparison: Individual + All Ensemble Methods" ))write_csv (comparison_all_models, "results_massive/complete_model_comparison.csv" )# ============================================================================ # VISUALIZATION 1: R² RANKING (Horizontal Bar Chart) # ============================================================================ :: cli_h2 ("Visualization 1: Model R² Ranking" )<- plot_ly (%>% arrange (R_squared),x = ~ R_squared,y = ~ Model,type = "bar" ,orientation = "h" ,marker = list (color = ~ R_squared,colorscale = "Viridis" ,showscale = TRUE ),text = ~ round (R_squared, 4 ),textposition = "auto" ,height = 900 ,width = 1500 %>% layout (title = "Model Performance: R² Score Ranking (Higher is Better)" ,xaxis = list (title = "R² Score" ),yaxis = list (title = "Model" ),margin = list (l= 300 , r= 50 , t= 80 , b= 50 , pad= 4 ),paper_bgcolor= 'white' # ============================================================================ # VISUALIZATION 2: Multi-Metric Comparison (Grouped Bar Chart) # ============================================================================ :: cli_h2 ("Visualization 2: Multi-Metric Comparison" )<- comparison_all_models %>% pivot_longer (cols = c (RMSE, MAE, R_squared), names_to = "Metric" , values_to = "Value" ) %>% mutate (Metric = factor (Metric, levels = c ("R_squared" , "MAE" , "RMSE" )))<- plot_ly (x = ~ Model,y = ~ Value,color = ~ Metric,type = "bar" ,text = ~ round (Value, 4 ),textposition = "auto" ,height = 900 ,width = 1500 %>% layout (title = "Model Performance: All Metrics (RMSE & MAE - Lower Better | R² - Higher Better)" ,xaxis = list (title = "Model" , tickangle = - 45 ),yaxis = list (title = "Value" ),barmode = "group" ,margin = list (t = 80 , l= 100 ),paper_bgcolor= 'white' # ============================================================================ # VISUALIZATION 3: SPEED vs ACCURACY TRADE-OFF (Scatter/Bubble Chart) # ============================================================================ :: cli_h2 ("Visualization 3: Speed vs Accuracy Trade-off" )<- data.frame (Model = comparison_all_models$ Model,InferenceSpeed = c (5 , 1 , 3 , 1.5 , 1.5 , 2 , 1.5 , 2.5 , 2.5 , 2.5 , 2.5 ) # Adjust to match models <- comparison_all_models %>% left_join (model_speed, by = "Model" )<- plot_ly (x = ~ InferenceSpeed,y = ~ R_squared,size = ~ MAE,color = ~ RMSE,text = ~ Model,mode = "markers" ,marker = list (sizemode = "diameter" ,sizeref = 2 * max (comparison_with_speed$ MAE) / 40 ^ 2 ,colorscale = "Viridis" ,showscale = TRUE height = 700 ,width = 1300 %>% layout (title = "Speed vs Accuracy Trade-off (Bubble size = MAE)" ,xaxis = list (title = "Inference Speed (seconds, lower = faster)" ),yaxis = list (title = "R² Score" ),coloraxis = list (colorbar = list (title = "RMSE" )),margin = list (l= 100 , r= 100 , b= 80 , t= 100 ),paper_bgcolor= 'white' # ============================================================================ # VISUALIZATION 4: USE CASE RECOMMENDATIONS (Heatmap) # ============================================================================ :: cli_h2 ("Visualization 4: Model Recommendations by Use Case" )<- data.frame (Model = comparison_all_models$ Model,"Fast Predictions" = c (5 , 4 , 2 , 4 , 4 , 3 , 4 , 3 , 3 , 3 , 3 ),"Highest Accuracy" = comparison_all_models$ R_squared * 10 ,"Robustness" = c (4 , 5 , 3 , 5 , 5 , 4 , 5 , 5 , 5 , 5 , 4 ),"Research" = c (3 , 4 , 4 , 4 , 5 , 5 , 4 , 5 , 5 , 5 , 5 ),"Real-time Apps" = c (5 , 4 , 2 , 4 , 4 , 3 , 4 , 3 , 3 , 3 , 3 ),"Regulatory" = c (3 , 4 , 2 , 4 , 5 , 4 , 4 , 4 , 4 , 4 , 4 )<- use_cases %>% pivot_longer (cols = - Model, names_to = "UseCase" , values_to = "Score" )<- plot_ly (x = ~ UseCase,y = ~ Model,z = ~ Score,type = "heatmap" ,colorscale = "RdYlGn" ,text = ~ round (Score, 2 ),texttemplate = "%{text}" ,height = 700 ,width = 900 %>% layout (title = "Model Suitability for Different Use Cases (1-10 scale)" ,xaxis = list (title = "Use Case" ),yaxis = list (title = "Model" ),margin = list (l= 150 , r= 50 , t= 80 , b= 80 ),paper_bgcolor = 'white' # ============================================================================ # VISUALIZATION 5: TOP 5 MODELS DETAILED COMPARISON # ============================================================================ :: cli_h2 ("Visualization 5: Top 5 Models - Detailed Metrics" )<- comparison_all_models %>% slice (1 : 5 ) %>% pivot_longer (cols = c (RMSE, MAE, R_squared),names_to = "Metric" , values_to = "Value" )<- plot_ly (x = ~ Model,y = ~ Value,color = ~ Metric,type = "scatter" ,mode = "lines+markers" ,height = 700 ,width = 1200 %>% layout (title = "Top 5 Models: Detailed Metrics Comparison" ,xaxis = list (title = "Model" ),yaxis = list (title = "Value" ),hovermode = "x unified" ,margin = list (l= 150 , r= 50 , b= 80 , t= 80 ),paper_bgcolor = 'white' # ============================================================================ # FINAL TEXT OUTPUTS # ============================================================================ <- comparison_all_models %>% slice (1 ):: cli_h2 ("🏆 Best Overall Model Performance" ):: cli_alert_success (paste0 ("Champion Model: " , best_overall$ Model)):: cli_ul (c (paste0 ("R² Score: " , round (best_overall$ R_squared, 4 )),paste0 ("RMSE: " , round (best_overall$ RMSE, 4 )),paste0 ("MAE: " , round (best_overall$ MAE, 4 ))# ============================================================================ # SAVE RESULTS # ============================================================================ write_csv (comparison_all_models, "results_massive/final_complete_model_comparison.csv" ):: cli_h1 ("✅ PRODUCTION PIPELINE COMPLETE WITH VISUALIZATIONS!" ):: cli_alert_success ("Generated 11 models total:" ):: cli_li ("3 Individual models (RF, XGB, DNN)" ):: cli_li ("4 XGBoost+RF Ensemble variants" ):: cli_li ("4 Multi-model Ensemble variants" ):: cli_alert_success ("5 Interactive Visualizations generated:" ):: cli_li ("1. R² Ranking Chart" ):: cli_li ("2. Multi-Metric Comparison" ):: cli_li ("3. Speed vs Accuracy Trade-off" ):: cli_li ("4. Use Case Recommendation Heatmap" ):: cli_li ("5. Top 5 Models Comparison" ):: cli_alert_info (paste0 ("🎉 Best Model: " , sub ("🏆 " , "" , best_overall$ Model), " (R² = " , round (best_overall$ R_squared, 4 ), ")" ))``` ```{} ``` `data_massive/` and `results_massive/` directories.